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A mathematical model of solidification, that simulates the formation of channel 
segregates or freckles, is presented. The model simulates the entire solidification process, 
starting with the initial melt to the solidified cast, and the resulting segregation is predicted. 
Emphasis is given to the initial transient, when the dendritic zone begins to develop and the 
conditions for the possible nucleation of channels are established. The mechanisms that lead 
to the creation and eventual growth or termination of channels, are explained in detail and 
illustrated by several numerical examples. 

A finite element model is used for the simulations. It uses a single system of 
equations to deal with the all-liquid region, the dendritic region, and the all-solid region. 

The dendritic region is treated as an anistropic porous medium. The algorithm uses the 
bilinear isoparametric element, with a penalty function approximation and a Petrov-Galeridn 
formulation. 

The major task was to develop the solidification model. In addition, this report 
briefly describes other tasks that were performed in conjunction with the modelling of 
dendritic solidification. 
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random number 

average total concentration of solute, 
heat capacity 

concentration of solute in the liquid 

interdendritic solute concentration from the liquidus line in the phase 
diagram 

local solute concentration in the solid • 

reference solute concentration (concentration in the liquid before solidification) 

primary dendrite arm spacing 

solute diffusivity in the liquid 

Darcy number in the x and z direction, respectively 

gravitational acceleration 

reference thermal gradient 

local mesh length in the direction of flow 

reference length 

z coordinate for the top of the container 

auxiliary variable for total solute concentration in the solid 

partition ratio 

permeability in the x and z direction, respectively 
latent heat 

first and second nondimensional latent heat, respectively 

slope of the liquidus line in the phase diagram 

unit outward normal vector 

shape function 

pressure 

hydrostatic pressure 

Petrov-Galerkin perturbation function 

Prandtl number 

prescribed boundary heat flux 

cooling rate prescribed at z = 0 

solutal and thermal Rayleigh number, respectively 

Schmidt number 

time 

temperature 
eutectic temperature 
initial temperature 
liquidus temperature 
reference temperature 
time step 

x component of the superficial velocity 
x component of the pore velocity 
reference velocity 
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v magnitude of liquid velocity 

solidification velocity 
z component of the superficial velocity 
z component of the pore velocity 
width of the container 

k Petrov-Galerkin weighting function 

z coordinates 

thermal diffusivity 
Petrov-Galerkin coefficients 

, /3 S thermal and solutal expansion coefficient, respectively, 

cell Reynolds number . 

7 3 cell Peclet number in the energy and 

solute concentration equation, respectively 
domain boundary 
convergence tolerance 

angle of the gravity vector with respect to the x axis 
penalty parameter 
kinematic viscosity 
density 

reference density 
reference time 

volume fraction of interdendritic liquid 
global domain 
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I. INTRODUCTION 


An increasing number of technological applications requires the operation of critical 
mechanical components under severe conditions of high temperatures and stresses. If 
improperly manufactured, these components are subject to creep fracture and thermal fatigue 
failures, which are almost always associated with grain boundaries that are transverse to the 
applied stress. The modem day directional solidification technique provides an effective 
means of controlling the grain shape, producing a columnar microstructure with all the grain 
boundaries running in the longitudinal direction of the casting. This greatly improves 
material performance at elevated temperatures. Still better properties can be obtained by 
casting a single crystal (i.e., a dendritic monocrystal), in which only one columnar grain is 
allowed to grow and form the main body of the casting [1]. 

Without proper control, directionally solidified castings are not free of defects. 
Particularly troublesome are localized segregates at the macroscopic scale, which are found 
in many directionally solidified (DS) alloys, including some superalloys. These segregates 
are observed as long and narrow trails, aligned parallel to the direction of gravity in DS 
castings, and are enriched in the normally segregating elements and depleted of the inversely 
segregating elements, i.e., their composition is shifted toward the eutectic composition. In 
horizontally solidified ingots of steel similar defects are known as "A" segregates, while in 
DS castings they have a more pronounced channel shape and are termed "freckles." This 
non-uniformity of composition is highly undesirable because the resulting variation in 
physical properties within the casting can lead to inferior performance of the components 
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manufactured from the casting. In ingot production, an excessive number of defects can 
require a large amount of cropping, at a considerable cost of energy and material. 

Many analytical and experimental works have been done in the recent years in order 
to observe, explain and prevent the formation of both "A" segregates [2-4] and freckles 
[5-13] in directionally solidified alloys. This paper addresses the latter type of defects, i.e . , 
freckles when an alloy is cooled from below. 

The opacity of metals prevents direct observation of the nucleation and growth of 
channels during solidification. Observations are usually done by quenching the ingot at 
different stages of solidification and analyzing the solidified macrostructure. Much of what 
is known about channel dynamics has been learned from the transparent analogue NH 4 C1- 
H 2 0 system [5,8]. Here, it is clearly seen that freckles are a direct consequence of upward 
flowing liquid jets that emanate from within the mushy zone. In the case of a binary alloy 
this requires that the solute be less dense than the solvent if it segregates normally or more 
dense if it segregates inversely. 

Because a density inversion also occurs in metallic systems where freckles are 
observed, it was proposed by Copley et al. [5], and further supported in later works 
[4,6,7,8,11,14], that buoyancy driven convection is responsible for channel segregation. 
Although the influence of buoyancy effects in segregation seems to be evident, the large 
differences in thermal conductivities, solid-liquid densities and, in particular, the fraction of 
liquid in the dendritic structure can lead to a convection pattern that is very different from 
that observed in the water-based mixtures, as it has been reported in experiments using 
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radioactive tracer techniques [14,15] and in an analytical model of thermosolutal convection 
in dendritic alloys [16]. 

Solidification simulation is motivated by the need to understand the basic' mechanisms 
of the segregation of alloy components. The earliest works considered only the solute 
conservation equation [17-19]. Soon thereafter, natural convection of the interdendritic 
liquid [20-22] and in the all-liquid region [23] were studied. Thermosolutal convection, 
however, was not considered, so early models were not capable of predicting thermosolutal 
instabilities at the solidifying interface that result in segregation defects in solidified alloys. 

Models that incorporate the effect of thermosolutal convection have been developed 
by other researchers as well, e.g., Bennon and Incropera [3,24-26], and Beckermann and 
Viskanta [27-28] simulated horizontal solidification of aqueous solutions of ammonium 
chloride. 

In our work, we have used a finite element algorithm to calculate macrosegregation 
and the formation of channels and freckles in Pb-Sn alloys. The algorithm has also been 
used to reproduce a calculation for horizontal solidification of an NH 4 C1-H 2 0 solution 
presented in Ref. [3]. Calculations presented herein show that the present model is capable 
of capturing the formation of channels and freckles during vertical solidification of alloys. 


II. SOLIDIFICATION MODEL 

We consider an alloy melt undergoing directional solidification under the following basic 
assumptions: 

1. Only solid and liquid phases may be present, i.e., no pores form. 
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2. The flow is two dimensional and laminar, and the solid phase is stationary. 

3. The solid and liquid phases have equal thermal properties and densities. 

4. There is no diffusion of solute in the solid phase. 

5. The thermal properties are constant, and the Boussinesq approximation is 
made; hence, the density is constant except in the body force term of the 
momentum equations. 

6. The dendritic region (often called the mushy zone) is treated as a porous 
medium with anisotropic permeability where the fluid velocities are the 
superficial velocity components. 

With these assumptions, the basic nondimensional equations of momentum, 
continuity, energy, and solute transport can be written [29] as 
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The variables are defined in the nomenclature. 

In the mushy zone, where 4> < 1, we use u and w as the superficial velocities, i. e . , 


u = <t>u p w = 4>w p (6) 

where u p and w p are the components of the pore velocity. 

The nondimensional numbers are defined by Pr = v/a, Da x = K x /H 2 , Da? - K z /H 2 , 

R t = 03 , r gGH 4 )/vot, R s = WsgC^H^/vD, Sc = v/D, and L = L/cGH. Note that His the 
reference length, G is a reference thermal gradient, and C M is the reference concentration. 
Implicit in Eqs. (1) and (2) is a constitutive relation of the form 

P = P 0 [1 ~ Pj (T - T 0 ) - 0 S (C - CJ] (7) 

in the body force terms, where p 0 and T 0 are reference density and temperature, respectively, 
and C^, the reference solute concentration, is the concentration in the liquid before 
solidification. The nondimensionalization given above also uses a reference velocity U = 
a/H, a reference time r = H 2 /v, and a reference pressure p 0 U 2 . The temperature and solute 
concentration in the liquid are nondimensionalized using (T - Tq)/GH and C/C*,, 
respectively, and the total solute concentration by C/C 
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The momentum equations above are similar to the ones used by Beckermann and 
Viskanta [27-28], but differ from the ones used by Bennon and Incropera [3,24-26]. Readers 
are referred to Ganesan and Poirier [30] and Nandapurkar et al. [31] for more details 
pertaining to the momentum equations. 

To complete the model, we require a relation between. the total solute concentration C 
and the solute concentration in the liquid C, of the form 

C = <f>C, + (1 - 4>)C S (8) 

where C s is the average local solute concentration in the solid. The model allows for 
microsegregation in the solid; hence, C s is not uniform but is given by 

l 

C s = _L [ kC, d4> (9) 

1_ * i 

where k is the equilibrium partition ratio defined as the ratio of the concentration of solute at 
the solid dendritic interface to the solute concentration in the interdendritic liquid. In this 
work, we assume k is constant. 

At a given temperature in the mushy zone, the composition of the interdendritic liquid 
is nearly uniform, and the local liquid-solid interface is very near equilibrium [32]; hence, 
the composition of the interdendritic liquid is given by the liquidus line in the phase diagram 
of the alloy. Consequently, 

C t = C l (T) (liquidus line in the phase diagram) 0 < 1 
C l = C 0 = 1 

and Eq. (8) can be used to calculate the local volume fraction of liquid, <£. 
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In this work, we present two-dimensional simulations in rectangular domains (Fig. 1). 
The boundary conditions associated with Eqs. (1) and (2) are no-slip at solid boundaries; 
i.e., u = w = 0 at solid boundaries that include completely solidified regions in the domain. 
For solidification in a very tall container, the normal stresses along the top boundary are 
zero,- viz . , 

$“+tL-0 2 Pr— = <t>p z = H t (11) 

dz dx d z 

If we want the upper boundary to be a free surface, we assume that the surface is 
undeformable. The boundary conditions are then 

^ = w = 0 z = H t (12) 

dz 

The boundary conditions for temperature along any of the walls can be of two types: 
a prescribed heat flux, 

-P 0 ca— = q (13) 

on 

(along the boundary), where d/dn denotes the normal derivative in the direction of the 
outward unit normal to the boundary and q is the prescribed heat flux along that boundary; 
or a prescribed temperature written as 

T = T. - rt (14) 

(along the boundary), where T t is the initial temperature and r is a prescribed cooling rate. 
Any combination of the above-cited conditions can be imposed along the boundaries of the 
domain. 
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No transfer of solute mass is allowed at solid boundaries or at an undeformable free 


surface (it is assumed at the top of the container), i.e., 

dC, nt ->, 

4> l = 0 ( 15 ) 

bn 

along those boundaries. If the assumption of a very long container is made, we must require 
a balance of diffusive and convective, transport along that boundary so that 

D— + w(C z - CJ = 0 (1.6) 

dx 

at the top boundary. 

III. CALCULATION STRATEGY 

We adopt a strategy in which the equations are solved sequentially and an iteration is 
performed to achieve convergence at each time step. First, we combine Eqs. (8) and (9) and 
obtain 

C = <t>q + I (17) 

where I is the integral in Eq. (9). 

To advance from time t n , at which all conditions are known, to time t n+J = t n + A t, 
the following steps are taken, where the dependent variables are computed using the latest 
available values of all other variables on which they depend. 

1. At time t n , parameters u n , w n , T n , etc., are known. 
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Advance to time step t n+1 = t n + At. Set i — 0, u ^ = u n , = w n , 

T^ x = T n , etc., where i is an iteration index and the subindex denotes the 


time level. 


Compute u l n * + \ and w' n l\ from Eqs. (l)-(3) and (7). 
Compute 7^ +1 from Eq. (4). 


From Eq. (10), set 



(13) 


Calculate <t>*\ 


from Eq. (17) in the form 
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< 1, compute I n+1 


(19) 


Compute C' n + + \ using Eq. (5). 

Recalculate <t>*\ using Eq. (19). 

If | Cm - Ci 1 < e (tolerance), then 

__ 1 + 1 1 + 1 rp __ rpi + l 

^n+1 “ ^n + l> "rt + 1 ^n+l> *h+1 


etc.', n - n + 1; go to step 2. 


Otherwise, set i = i + 1 and go to step 4. (Here, || / || denotes the 
Euclidean norm.) 



Steps 4-10 are repeated iteratively during each time step until convergence of 4> is achieved. 
The velocities are calculated only once per time step because they have negligible sensitivity 
to small changes in the rest of the dependent variables [29]. At the end of each time step, 
the average concentration in the solid C s is calculated from Eq. (9). 

A. Remelting 

The integral 7 in Eq. (17) is computed incrementally by adding the increment 
corresponding to a change in <jj to the previous value of 7. The increment can be positive or 
negative, depending on whether the material undergoes solidification or remelting. More 
details can be found in Ref. [33] . If solidification has occurred and the increment is positive, 
7 is calculated directly. If there is remelting, however, 7 must be obtained from the 
solidification history, with values of 4> and 7 saved from the previous time steps. To 
alleviate the excessive amount of storage required, <f> and 7 are not saved at every time step. 
Instead, 7 is stored at increments of A<£ = 0.01, and linear interpolation is used for other 
values of <£. 

B. Energy Equation 

A simple analysis shows that the algorithm for calculating the temperature and the 
volume fraction liquid is only conditionally stable and, unfortunately, stability only holds for 
values of latent heat that are much smaller than in metallic alloys. In order to obtain a stable 
algorithm, the energy equation, Eq. (4), was reformulated to make the latent heat term 
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implicit. To eliminate <f> from Eq. (4), we first differentiate Eq. (17) with respect to time, 
yielding 


8C = 8C, 

dt * dt 


+ c t — 

1 dt 



( 20 ) 


• If we approximate the liquidus line by a straight line with slope m (m < 0), the 
temperature is related to the solute concentration in the interdendritic liquid by 


T = T m + jftC, (21) 

where T M is the melting point of the solvent and m = itiCq/GH. Equations (5), (20), and 
(21) can then be combined with Eq. (4) to give 
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( 22 ) 


where L is the second nondimensional latent heat defined as L = L/mcC 0 . Note that when 
remelting occurs, we cannot make the assumption that d/dr[ j kC t d(f> ] = -kC l d<f>/dt, as 


was done, for example, in Ref. [17]. An explanation of this point is given by Felicelli [29] 
and by Felicelli et al. [33]. 

The fact that Eq. (22) rather than Eq. (4) must be used to calculate the temperature 
introduces a computational inconvenience because Eq. (22) is nonlinear and must be 
reconstructed and solved again at every time step, increasing the CPU time significantly. It 
must also be pointed out that Eq. (21) is not uniformly valid throughout the domain because 
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it does not reduce to Eq. (4) when <f> = 1. Therefore, the terms containing L apply only to 
the elements in the mushy zone, where C ; = C L (T). 

C. Solidification at Eutectic Temperature 

Binary alloys solidify over a range of temperatures, and the temperature is governed by 
the energy equation, Eq. (4). If we look back at the computational strategy, the algorithm 
proceeds as follows: 

1. A decrease in the temperature of a volume element in the mushy zone 
automatically results in an increase in the solute concentration in the 
interdendritic liquid, which is given by Eq. (18). 

2. An increase in the concentration in the interdendritic liquid produces a 
decrease in the volume fraction of liquid according to Eq. (19). 

This means that the temperature is the variable that drives the solidification process. When a 
volume element reaches the eutectic temperature, however, the remaining fraction of 
interdendritic liquid solidifies at constant temperature, assuming no undercooling of the 
eutectic reaction. Consequently, at this point the solidification algorithm must be modified. 
In formulations based on enthalpy, such as the one used by Bennon and Incropera [24-26], 
the enthalpy is still valid at the eutectic temperature. 

In some ways, enthalpy formulations are convenient, but they require a transformation 
to obtain the temperature. Unfortunately, an algebraic transformation relating enthalpy to 
temperature is only possible for simple thermophysical descriptions, such as constant specific 
heat and complete diffusion in the solid. The assumption of complete diffusion in the solid 
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yields a direct relation between <f> and T that, when combined with the additional assumption 
of linear dependence of the phase enthalpies with temperature, yields an algebraic equation 
for the temperature as a function of the mixture’s enthalpy. This property was Exploited by 
Bennon and Incropera [3], who used an enthalpy-based energy equation. 

In order to retain the temperature formulation, an alternative is followed to solve for 
solidification at constant temperature. When a node in the domain reaches the eutectic 
temperature, the energy equation, Eq. (4), is used to calculate the volume fraction of liquid 
while setting the time derivative of T equal to zero to effect solidification at constant 
temperature. Equation (4) becomes 


30 = _1 d 2 T ^ B 2 T _ u dT _ w dT 

dt L fix 2 dz 2 U dx W dz 


(23) 


Equation (23) is solved only at the points solidifying at constant temperature; the complete 
energy equation is solved elsewhere. Physically, the eutectic isotherm advances smoothly 
through the domain; that is, a coordinate does not stay at the eutectic temperature for a finite 
period of time. The present model, however, treats volume-averaged quantities in a porous 
medium, so a volume element or nodal point stays at constant temperature until it is 
completely solidified. After that, its temperature is calculated with Eq. (4), which reduces to 
the conduction equation in the solid. 


D. Time Step and Mesh Size 

If thermosolutal effects are important, the relevant length scale to be resolved is D/V, 
where Vis a characteristic solidification speed. For the system under consideration, i.e . , 
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slowly solidifying Pb-Sn alloys, this length scale is typically of the order of 300 pm and 
requires a very fine mesh spacing in the vertical direction. In the horizontal direction, 
channels in the mushy zone can be narrow, comparable to the order of D/V. Finally, we 
work with uniform meshes because the location at which channels form cannot be predicted 
a priori; as a result, calculations are limited to the rather small domains. 

A second consideration for selecting the spatial resolution is that there be enough nodes 
to adequately resolve the distribution of the volume fraction of liquid $ in the mushy zone. 
This consideration is discussed later in relation to the mesh spacings used in the calculated 
results. 

The choice of the time step is more complicated because of the nonlinear nature of the 
problem and the need to satisfy several stability constraints in the calculations. So far, the 
calculations are performed using a constant time step; hence, the most stringent of the 
following stability criteria must be satisfied. 

1. The convective terms are treated explicitly in Eqs. (1), (2), and (4), so the 
stability condition of the form At { < ( \ u \ /Ax + \ w | /A z)' 1 applies 
everywhere in the computational mesh. 

2. The explicit treatment of the body force terms gives A t 2 ^ (R T + R s /Sc)''. 

3. A third consideration can be added, namely, that it is more efficient to choose 
a time step for which no more than three iterations are needed for convergence 
within the time step. 

The above criteria are not optimal bounds. However, choosing At < min (A t u At 2 ) 
gives a good estimate. An added difficulty is that A t x can vary significantly as the 
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solidification proceeds. Hence careful monitoring of the flow field is required to avoid an 
overly conservative time step. 


IV. FINITE ELEMENT METHOD 

The computational algorithm makes use of a standard penalty function formulation 
because the pressure is not needed. The pressure in Eqs. (1) and (2) is eliminated using the 
pseudo-constitutive relation 



(24) 


where p s is the hydrostatic pressure and X is a large "penalty" parameter. This method has 
been discussed in detail by Marshall et al. [34], Heinrich and Marshall [35], and Heinrich 
and Yu [36] in the context of buoyancy-driven flows. Theoretical aspects have been studied 
by Oden [37]. Equation (24) is substituted into Eqs. (1) and (2), and the continuity equation, 
Eq. (3), is eliminated. As a result, the pressure is not calculated but, if needed, it can be 
recovered a posteriori by solving a Poisson equation (see [38, 39]). 

The present scheme is based on the bilinear Lagrangain isoparametric element. The 
convective terms are dealt with using a Petrov-Galerkin formulation in which the weighting 
function is perturbed in the convective term. If we denote the shape function corresponding 
to node i by N t , the perturbed weighting function takes the form 


W* = N t + a K P t K= 1,2,3, < 25 > 

where K = 1 corresponds to Eqs. (1) and (2), K = 2 to Eq. (22), and K = 3 to Eq. (5). 
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The parameters a K are given by 


where 


a K ~ co ^h y k — 
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(26) 
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is the solute concentration cell Peclet number, v = | \ju 2 + w 2 | is the local fluid speed, 
and h is the local element length in the direction of flow. Details on calculating h are given 
in Ref. [36]. 

Denoting the domain by fl and its boundary by T, the weak forms of Eqs. (1) and (2) 
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where the left-hand sides are treated implicitly and the right-hand sides are evaluated 
explicitly with the latest values of the dependent variables. All variables are interpolated 
using the bilinear isoparametric shape functions; i.e., for a generic function f(x, z) we write 


fix, z) = £ Nfic, z) (32) 

I 

where i ranges over the number of nodes in the mesh and f t = fix it z,) is the value of the 
function at note /. 

The weak form of the energy equation, Eq. (22), is 
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In Eq. (33), the line integral is taken only over those portions of the boundary where a flux 
condition as given by Eq. (13) is prescribed. 

The weak form for the solute concentration equation is given by 
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Sc bn 

where the line integral is only calculated along the top boundary when the long-container 
assumption is made and Eq. (16) is used. In the all-liquid region, C = C h and the left-hand 
side is treated implicitly. In the mushy zone where <f> < 1, however, we have C ^ C,, and 
the second term on the left-hand side is calculated explicitly using the latest values of C,. 

After interpolating the functions using Eq. (32), all integrals can be evaluated, and 
three sets of linear algebraic equations result: 


M iV + K t v = - Qj(v, <(>, T, C) 


(35) 


M QT * K 2 (4>, C,)T = Q 2 (v, 0, T, Q 


(36) 


M 3 C + K 3 (4>)C, = - Q 3 (v, <}>, C,) 


(37) 


In the equations above, v = (u l w l u 2 w 2 . . . u N w N ) T , where N is the number of nodes and u t 
and w,- are the velocity components at the nodes; also, = (<^> 1 0 2 • • • 4 > n ) t , T = 

(T x T 2 . . . T n ) t , C x = (CiCi 2 . . . C,/, and C = (CjC 2 . . . C N ) T . The mass matrices 
and M 2 are lumped diagonal matrices (see Refs. [36, 40, and 41] and are given by 


18 



Klj = }„ Pf w, 
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(39) 


(i = 1, 2, . . . , AO 

The mass matrix M 3 in the solute concentration is treated consistently and is given by 

[m 3 lj = | fl Pr N^jcm i, j = 1, . . . , N 
The stiffness matrices K ; are given by 
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The vectors Q ; are obtained by replacing Eq. (32): 
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The time integration of the systems of ordinary differential equations represented by 
Eqs. (35)-(37) is performed using a backward implicit scheme. Specifically, discretize the 
time derivative of a scaler function /using 


/ ^ 


= r +l -r 


At 
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V. NUMERICAL RESULTS 


Calculations were performed for a Pb-10 wt pet Sn alloy under a variety of cooling 
and boundary conditions. The computational domain is shown in Fig. 1 and consists of a 
rectangular region of width W and height H T , discretized with a uniform mesh of 20 elements 
in the x (horizontal) direction and 30 elements in the z (vertical) direction. This mesh was 
selected after several preliminary calculations were performed with different resolutions to 
assess the sensitivity of the results to the mesh size. The chosen mesh was found to capture 
all the dynamic features of the flow at a reasonable computation cost, whereas coarser 
meshes were unable to reveal the formation of the channels. A uniform resolution was used 
because the locations of channels and regions of macrosegregation were found to be 
unpredictable. 

In all calculations, the channels turned out to be very narrow, comparable to one to 
several primary dendrite arm spacings (which is also comparable to the solute distribution 
decay scale D/V). The mesh size has to be kept comparable to this scale if proper resolution 
of the channels is desired. This requirement imposed a limitation on the size of the domain, 
in order not to make a single computation excessively expensive. Fortunately, all the effects 
associated with channel creation and growth could be observed in calculations even though 
the domains were rather small. Furthermore, when different sizes were used, it was found 
that channels developed under the same thermal conditions irrespective of size of domain. 

A reference length scale H = 300 fim was chosen for all the calculations. The 
following nondimensional parameters were also used: 
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Ra T = 3.4 x 10' 3 Ra* = 1.86 x 10 4 

Pr = 2.3 x 10; 2 Sc' = 8.23 x 10 1 

L/tGH = 7.5 X 10 2 

The reference time, r, was 0.65 s, and the reference velocity, U, was 27 mm-s 1 . The 
reference temperature was taken to be the freezing point of the Pb-10 wt pet Sn alloy 
(T 0 = 577 K), and the reference concentration was C 0 = 10 wt pet Sn. 

The first calculation considers a container of width W = 10 mm with a domain height 
of H t = 20 mm. Zero tangential stress and zero vertical velocity are imposed on the top 
boundary (z = H T ), and no-slip at the other three boundaries. The alloy is initially all liquid 
of concentration C 0 = 10 wt pet Sn, and has a linear temperature distribution varying from 
T 0 (the melting point) at the bottom boundary to T 0 + GH T at the top boundary, with G = 
1000 K-m' 1 . The side walls are insulated, and a vertical gradient of dT/dz - G is imposed 
at the top boundary. A time dependent boundary condition is prescribed at the bottom of the 
form: 

T(x,Q,t) = T 0 - rt 

where r is the cooling rate (0.0783 K s' 1 ). With this value of r, the initial solidification rate 
was approximately 0.042 mm s' 1 . The values of the cooling rate, r, and the thermal 
gradient, G, were selected from the thermal history of ingots of Pb-10 Sn alloy that exhibited 
channel segregates in the experiments of Sarazin and Hellawell [7). 

The thermodynamic and transport properties used in the calculations are given in 
Table I. The functions used for the anisotropic permeability are given in the Appendix. 
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Table I. 


Thermodynamic and Transport Properties Used for Calculations 
in Pb-10 wt pet Sn. Taken from Ref. [33]. 


Property 


Reference concentration (C 0 ), 10 wt pet Sn 

Reference temperature (Tq), 577 K 

Equilibrium partition ratio (k), 0.31 

Melting point of lead (T M ), 600 K 

Slope of liquidus (m), -2.33 K (wt pet)' 1 

Density [p 0 = p(C 0 ,T 0 )], 1.01 X 10 4 kg-m' 3 

Thermal expansion coefficient (/3 r ), 1.2 X 10' 4 K* 1 

Solutal expansion coefficient (/3 C ), 5.15 x 10' 3 (wt pet) 1 

Kinematic viscosity (j> 0 ), 2.47 x 10' 7 m 2 s 1 

Latent heat (L); kJ kg' 1 : 

At (Cq,T q ), 26 
At (C e ,T e ), 56 
Used in calculation, 37.6 
Heat capacity (£); kJ kg' 1 K' 1 : 

Liquid at (C 0 ,r 0 ), 0.161 
Solid at (C e ,T^), 0.173 
Used in calculation, 0.167 
Thermal conductivity (k); kW m' 1 K' 1 : 

Liquid at C 0 ,T 0 ), 0.0167 
Mixture at ( C E ,T 0.0198 
Used in calculation, 0.0182 
Thermal diffusivity (a = K/p 0 t), 1.1 x 10' 5 m 2 s' 1 
Solutal diffusivity (D), 3 x 10" 9 m 2 s' 1 


Figure 2 shows calculated results after 3 min. The contours of equal fraction liquid 
(Fig. 2a) show a very narrow region of high fraction liquid at the center of the casting, 
indicating the presence of a freckle or a channel. Note that the mesh size is not sufficiently 
small to resolve an all-liquid jet. Nevertheless, it can be seen from the streamlines (Fig. 2b) 
that the flow is upward within the channel and extends into the bulk liquid as a column of 
rising fluid. Although the streamlines penetrate into the mushy zone, the strength of the 
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convection rapidly diminishes within the mushy zone where the permeability is low. 

Figure 2c shows that the fluid within the channel is richer in solute than the surrounding 
fluid. 

A more detailed view of a freckle is obtained in the next example (Fig. 3), where a 
very small container of dimensions 2.5 mm by 4.5 mm was used, together with a slower 
initial solidification rate of V = 0.0139 mm-s* 1 . The width of the channel is approximately 
equal to the reference length scale (H = 0.3 mm). Note the steep gradient of fraction liquid 
on the channel walls, indicating that the channel is surrounded by a region of low volume 
fraction liquid. The flow is upward inside the channel, and two cells have developed next to 
the mold walls (Fig. 3b), where there is also upward flow and a higher than normal fraction 
of liquid. 

In the examples of Figs. 2 and 3 the channels in the centers were induced by 
introducing an initial perturbation in the concentration of the melt along the vertical 
centerline and letting the system evolve thereafter. The form of the perturbation was 

C, = 1.01 C 0 for x = W/2. 

In this way, an initial upward flow at the center of the container is established. The situation 
is similar to experiments done by Sample and Hellawell [8] in NH 4 C1-H 2 0, in which they 
created channels in the mushy zone by inducing an upward flow in the liquid just ahead of 
the dendritic tips by suction up through a capillary tube. 

In the next calculation no perturbation is introduced, and the system is left to evolve 
undisturbed. The width of the mold is 5 mm, and the height of the computational domain is 
H t = 10 mm. No-stress conditions are imposed on the top boundary, simulating an 
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infinitely high column of liquid above the mushy zone. The cooling rate is r = 0.0167 K 
s' 1 , corresponding to an initial solidification rate of V = 0.0139 mm s" 1 (50 mm h' 1 ). All 
other boundary conditions were the same as in the examples of Figs. 2 and 3. The results of 
this calculation are given in Figs. 4 and 5. 

After 2 minutes, two channels begin to form (Fig. 4a) at locations corresponding to 
zones of upward flow of four convective cells (Fig. 5a). The channels, however, do not 
keep growing upward, but they turn toward the sides and continue to grow along the walls of 
the mold (Figs. 4b through d), following the path of upward convective flow (Figs. 5b 
through d). An additional channel can be seen in Figs. 4c and 5c that appears as a pocket in 
the upper part of the mushy zone (indicated by arrows). However, the channel rapidly closes 
leaving a pocket of liquid surrounded by a dendritic network with a greater volume fraction 
of solid (Fig. 4d). Abrupt transitions in the volume fraction of liquid surrounding the pocket 
or the channels are indicated by the dark regions of closely packed contour lines. 

A plot of total concentration of the partially solidified casting (Fig. 6) shows that 
within the mushy zone there is obvious positive macrosegregation in the channels along the 
wall and within the now isolated pocket. 

The effect of heat flow between the casting and the furnace is studied in the next 
example. The two calculations are shown in Figs. 7a and b for a domain with a width of 
5 mm and a no-stress top boundary at 10 mm as in the previous example. The cooling rate, 
however, is only 8.33 x 10' 3 K s' 1 . In Fig. 7a heat is also extracted from the casting by 
imposing a temperature gradient | dT/dx \ = 0.1 K mm' 1 at the side walls. Except for 
channels at the walls, the mushy zone adopts a concave shape toward the bulk liquid. Only 
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channels at the walls develop and no interior channels or pockets are observed. In Fig. 7b, 
heat is added to the casting by imposing | dT/dx \ = 0.3 K mm' 1 at the side walls and by 
specifying a lower, temperature at the center of the base (increasing linearly toward the 
sides). Apart from the channels, the overall shape of the growth front is convex to the 
liquid; now channels develop in the interior of the container, although the channels at the 
walls still remain. 

It is evident from the above example that channels tend to develop in the leading part 
of the growth front. This fact has been observed in experiments with NH 4 C1-H 2 0 mixtures 
[5,8]. Copley et al. [5] varied the configuration of their bottom chill so as to make the 
growth front either concave or convex to the liquid. When it was concave, channels in the 
NH 4 C1-H 2 0 system tended to form on the outside, and when it was convex the channels were 
prevalent in the center. Similar results were obtained by Sample and Hellawell [8], who 
tilted the ingots. 

Based on the calculations presented, herein, it is also evident that there is a strong 
tendency for freckles to form at the walls of the mold. In order to investigate why a wall 
might be an attractive place for freckling, a numerical experiment was done using the same 
domain with insulated walls and prescribing a zero horizontal velocity along the vertical 
center line (/. e . , * = W/2). It is observed in Fig. 8a that a freckle forms at the center of the 
casting and keeps growing, supported by a column of upward flowing liquid, in a similar 
way as the channels along the walls (Fig. 8b). Two additional channels develop at positions 
intermediate between the center and the walls, but they do not persist (Fig. 8c) because they 
face an unfavorable flow pattern in the liquid zone. It is important to remark that in this 
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example the upward flow at the center of the container is not induced by an initial 
perturbation, but it arises naturally as a consequence of restricting the horizontal movement 
of fluid. A similar result was obtained by prescribing a highly anisotropic permeability 
(K x = 10' 3 K z ) along the vertical center line in the mushy zone, showing that a restriction to 
horizontal convection in the mushy zone can sustain a column of upward flow in the liquid. 

VI. DISCUSSION 

The results presented herein illustrate some of the most common effects observed in 
experiments on freckles. Several remarks can be made that are based upon the results. 

1 . The convection starts immediately, as soon as the liquid next to the base of the 
container begins to solidify. The convection cells nucleate near the advancing 
solidification front (/.£., at the tips of the dendrites). From the calculated 
temperatures and concentrations in this region, it is found that the vertical density 
gradient is positive there, i.e., it is gravitationally unstable. 

2. As the solidification proceeds, the convective cells do not keep organized. The 
flow pattern changes continually, driven almost entirely by the solute 
concentration field. The isotherms are practically flat because of the low Prandtl 
number. The convection in the bulk liquid penetrates deeply into the mushy 
zone, although the velocities decrease by 2 or 3 orders of magnitude below the 
upper 20 pet of the mushy zone. 

3. The first channels begin to form in regions of upward flow, between two 
convective cells ( e.g ., Fig. 2a). Because the solute concentration in the liquid is 
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higher at the bottom of the casting, the upward flow transports solute-rich liquid 
to the tips of the dendrites, decreasing the freezing point. Therefore, the growth 
front advances more slowly in regions where there is upward flow, it distorts, 
and a channel begins to develop backward from the tips toward the bottom of the 
mushy zone. This process is helped further by remelting of already existing solid 
in the mushy zone. The diffusion of solute is too slow to spread away from the 
column of rising solute-rich liquid, so remelting occurs in order to decrease the 
solute concentration in the liquid back to an equilibrium value that corresponds to 
the local temperature. In other words, without remelting a volume of rising 
interdendritic liquid is unable to maintain its thermodynamic equilibrium because 
heat diffuses much faster than the solute; remelting counteracts this by diluting 
the liquid with respect to the solute. 

4. The larger concentration of solute in a column of rising liquid tends to keep it 
rising because of buoyancy and double diffusion effects (assuming that the solute 
is lighter than the solvent, e.g. Sn in Pb-rich alloys). This keeps the channel 
growing vertically, until the pattern of convection changes, and the column of 
liquid emerging from the channel is disturbed. The flow change is caused by 
accumulation of solute in a nearby place, creating a competing channel. The old 
channel deviates from its vertical trajectory, induced by the upward flow at the 
new position, and eventually merges with the new channel or closes before 
reaching it, forming a pocket or streak. 
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5. In most of our calculations, we observe a strong accumulation of solute at the 
walls of the mold, leading to- long, well defined channels that dominate the flow 
and prevent other similar channels from growing in the interior of the casting. 
Channels that develop in the interior of the casting, either naturally or induced by 
a perturbation, are eventually absorbed by the channels at the walls. This 
predominance of freckles at the surface of the mold has been observed in 
Ni-based alloys [6]. The solute accumulates next to the walls because of the 
limitation of lateral transport; i. e . , it is harder to remove solute away from the 
wall because the horizontal component of the velocity is small. For high 
solidification rates (—4.2 X 10' 2 mm s' 1 ) the formation of channels at the walls 
can be delayed or prevented, because not enough time is allowed for solute 
accumulation. 

6. Long, vertical channels, similar to the ones at the surface of the casting, can be 
obtained in the interior if a restriction to lateral convection is introduced (Fig. 8). 
The purpose is to provoke an accumulation of solute that can sustain the channel 
growth. Such an accumulation could occur in regions of high anisotropy in the 
mushy zone (K x « K z ), as might be the case in grain boundaries, dendrite 
misalignments or other defects in the dendritic network of the mushy zone. 
Another mechanism that can favor the formation of long interior channels is the 
deposit of dendrite debris, produced by remelting or erosion of the dendrite 
arms, along the channel sides. This phenomenon has been observed in the 
NH 4 C1-H 2 0 system [5]. The debris grow as small randomly oriented grains that 
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entangle with the sides of the channel and reduce the horizontal permeability. In 
effect, the solute is confined to the channel. Because this decreases the lateral 
feeding of liquid, the channel activity is significant only in the neighborhood of 
the dendrite tips, and the liquid in the bottom part of the channel is essentially 
stagnant. 

VII. CONCLUSIONS 

A mathematical model of solidification of dendritic alloys with thermosolutal 
convection has been presented. The solidification is initiated from an all-liquid state, and the 
dendritic zone is allowed to grow as the volume fraction of liquid in the mushy zone adjusts 
according to local thermodynamic equilibrium conditions. Calculations were performed for 
the vertical solidification of a Pb-10 wt pet Sn alloy, with the following results: 

1. Channels form during the growth of the mushy zone. 

2. The shape of the channels vary from small pockets or short streaks of 
interdendritic liquid to long vertical penetrations. 

3. The formation of a channel occurs near the leading part of the growth front. 

4. The liquid within the channels is enriched in solute. 

5. The channels form in regions of upward flow. The mushy zone grows around 
the columns of rising liquid, forming the channels backward from the tips of the 
dendrites toward the bottom of the mushy zone. Remelting also takes place. 

6. The convection is driven mainly by the solute field.- Zones of solute 
accumulation are potential starters of channels. 
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7. Channels can be induced by establishing a column of upward flow in the bulk 
liquid, by introducing a restriction to the horizontal convection in the mushy 
zone, or by enhancing the vertical convection. 

8. There is a strong preference for channels to form at the mold walls. 

VIII. OTHER ACCOMPLISHMENTS 

The major emphasis in the research program was to model thermosolutal convection 
during the dendritic solidification of binary alloys, and the results of this effort are given in 
the preceding chapters. In this chapter, we summarize other accomplishments pertaining to 
the subject research. Details are not given here, but publications resulting from the work 
may be consulted. 

The two most important publications on the model, itself, are given as References 33 
and 46. The former appears in Metallurgical Transactions, which is widely read by the 
materials community, and the latter appears in Numerical Heat Transfer, where the 
mathematical details pertaining to the finite element algorithms and computational strategy 
peculiar to the solidification scenario are given. Our results were also presented at several 
conferences, all of which resulted in publications (References 47-50). 

In setting up the conservation equations for the mushy zone, we found it necessary to 
derive the momentum equation [30] and to numerically test our calculations [31] and contrast 
them to results calculated from a form of the momentum equation used by other researchers. 
This was necessary because in analyzing dendritic solidification * the mushy zone is modeled 
as a porous medium with a spatially varying fraction of liquid, unlike the uniform porosity in 
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most porous media. There were also fundamental issues, which had to be addressed, in both 
the energy equation and the solute conservation equation, resulting in the publication of 
Reference 51. 


APPENDIX- 

Functions Used for Permeability 

We assume directional solidification with dendritic columnar grains. The 
permeability is anisotropic with components K x and K y for flow perpendicular to the 
columnar dendrites and flow parallel to the columnar dendrites, respectively. 

In this work thermal conditions, corresponding to ingots studied by Sarazin and 
Hellawell [7], were used in the simulations, so we also used their reported primary dendrite 
arm spacing of 300 /xm. Lacking a value for secondary dendrite arm spacing, Eq. [30] from 
Poirier [42] was selected for K x ; it is 

K x = 7.08 X 10' 16 d 2 x °* <£ 3 - 32 t A1 ] 

with K x in m 2 and d x in /xm. 

Equation [Al] is based on a regression analysis of empirical data that cover the scope 
0.19 < <f) < 0.66 and 144 < d x < 420 /xm. 

For the regime with <f> > 0.66, we use 


d] 


Uq + 


0 

1 - 4 > 


1/4 


where % = -5.955449 x 10' 2 and a x = 5.652925 x 10" 2 . 


[A2] 
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Equation [A2] and the coefficients were determined by a linear regression of results 
given in Tables 1 and 2 of Sangani and Acrivos [43] for 0.6 < < 0.95, and calculated 

with their Eqs. (17) and (25) for 0.95 < $ < 0.99. Sangani and Acrivos [43] studied the 
slow flow past square arrays and triangular arrays of circular cylinders; Eq. [A2] is a 
compromise between the square array and the triangular array. 

Flow Parallel to Columnar Dendrites 

The empirical data are listed in Reference 42; the scope is 0.17 < 4> ^ 0.61 and 
28 < d x < 420 (im for which- the following empirical relationship was determined: 

— - = 3.75 x 10 -4 (f) 2 [A3] 

d\ 

There are no empirical data for cf> > 0.61, so we use the results of Drummond and Tahir . 
[44] and Ganesan et al. [45]. 

Drummond and Tahir [44] derived equations for laminar flow parallel to regular 
arrays of circular cylinders. Based upon values of drag forces given in their Table 5 and 
their equations for small values of (1 - <f>), permeabilities were calculated. 

Ganesan et al. [45] calculated the permeability for flows parallel to primary dendrites 
based upon actual microstructures. Their results agree very well with the results derived 
from Drummond and Tahir for g L > 0.55. Therefore, when <f> > 0.61 we use an equation 
given by them; it is 
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a 


^ = 

d\ 


In (1 - 4 >)~ 1 - b + 2(1 - <j>) - ij— 


[A4] 


where a and Z> are constants that depend on the array of the cylinders. In this work, we 
selected average values for triangular and square arrays: a = 0.07425 b = 1.487. 
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Fig. 1— Domain and coordinate system for vertical solidification. 
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BY 0.25 


Fig. 6-Total concentration (weight percent in solid plus liquid) showing macrosegregation in 
channels, corresponding to Figs. 4d and 5d (7 min). 
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Fig. 7— Channel locations by volume fraction liquid in a (a) concave (5 min) and (b) convex (20 min) mushy zone. 
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Fig. 8-Channel at the center of the casting formed by inserting a restriction on horizontal convection: (a) contours of fraction 
liquid at 10 min, (b) stream functions at 10 min, and (c) contours of fraction liquid at 15 min. 



